MODULE GMM_SE
!this module calculate the Standard Error of the GMM estimator
! decrease the theta by certain percentage
  USE GLOBVAR
  USE NRUTIL, only : ARRAY2D, ARRAY3D
  use comf_trans, only : real2ab_exp, ab2real_exp
  use MATRIX, only : Matrix_Inverse
  use DISTANCE

  IMPLICIT NONE
  
  character(4) :: gmmse_filename
  TYPE(ARRAY2D) :: gmmse_jacobian, gmmse_W, gmmse_W2, gmmse_XI2_PR
  TYPE(ARRAY3D) :: gmmse_simmoments
      
CONTAINS

  subroutine sub_GMM_SE(theta_)  !called in the master
    implicit none
    REAL(8), INTENT(IN) :: theta_(:)
    real(8) :: rprime(N_PAR)
    integer(4) :: i
  
    gmmse_simmoments%d1 = max(N_data_10, N_data_10_ss)
    gmmse_simmoments%d2 = MSM_Moments_Data%d2  !6 (no health) or 9 (with health) or 7 (pt only) or 10 (pt and health)
    gmmse_simmoments%d3 = N_PAR + 1 !the last one is the original one
    allocate(gmmse_simmoments%a(gmmse_simmoments%d1, gmmse_simmoments%d2, gmmse_simmoments%d3))
    gmmse_simmoments%a = 0.0d0
    
    ! moments:
    ! 1 lfpr_mean, 2 lnw_mean, 3 lnw_sd, 4 lnw_fe_mean, 
    ! 5 consumption  6 ss,  7 diff_lfpr_E2G, 8 diff_lfpr_G2B, 9 diff_lfpr_B2D, 7/10 pt
    ! 1.1 E to U, 2.1 U to E, 3.1 wage depreciation after unemp
    
    !keep original moments. This has to run immediately after the output
    gmmse_simmoments%a(1:gmmse_simmoments%d1,1:MSM_Moments_Sim%d2,gmmse_simmoments%d3)  &
          & = MSM_Moments_Sim%a(1:gmmse_simmoments%d1,1:MSM_Moments_Sim%d2)
     
    gmmse_W%d1 = MSM_wgt%d1        !#moments: N_moments_all
    gmmse_W%d2 = gmmse_W%d1        !#moments: N_moments_all
    allocate(gmmse_W%a(gmmse_W%d1, gmmse_W%d2))
    gmmse_W%a = MSM_wgt%a    !use the diagonal re-scaled weighting matrix
    
    ! OPTIMAL WEIGHT: only inverse of the diagonal.
    gmmse_W2%d1 = MSM_wgt%d1        !#moments: N_moments_all
    gmmse_W2%d2 = gmmse_W%d1        !#moments: N_moments_all
    allocate(gmmse_W2%a(gmmse_W2%d1, gmmse_W2%d2))
    gmmse_W2%a = 0.0d0
    do i = 1, gmmse_W2%d1
       if (gmmse_S%a(i,i) .NE. 0.0d0) gmmse_W2%a(i,i) = 1.0d0 / gmmse_S%a(i,i)  
    enddo
    
    gmmse_jacobian%d1 = N_moments_all     ! which is same as MSM_wgt%d1
    gmmse_jacobian%d2 = N_PAR
    
    allocate(gmmse_jacobian%a(gmmse_jacobian%d1,gmmse_jacobian%d2))
        
    gmmse_XI2_PR%d1 = gmmse_W%d1          !#moments: N_moments_all
    gmmse_XI2_PR%d2 = gmmse_XI2_PR%d1     !#moments: N_moments_all
    allocate(gmmse_XI2_PR%a(gmmse_XI2_PR%d1, gmmse_XI2_PR%d2))
    

  if (1>0) then  
    rprime = 1.0d0-1.0d-1
    write(gmmse_filename, "(A3)") '_N1'
    call sub_GMM_SE_get(theta_, rprime)
      
    rprime = 1.0d0-5.0d-2
    write(gmmse_filename, "(A4)") '_5N2'
    call sub_GMM_SE_get(theta_, rprime)
      
    rprime = 1.0d0-1.0d-2
    write(gmmse_filename, "(A3)") '_N2'
    call sub_GMM_SE_get(theta_, rprime)
  endif
        
    !deallocate variables
    deallocate(gmmse_simmoments%a)
    deallocate(gmmse_jacobian%a)
    deallocate(gmmse_W%a)
    deallocate(gmmse_W2%a)
    deallocate(gmmse_XI2_PR%a)
    
  end subroutine sub_GMM_SE
  
  
  subroutine sub_GMM_SE_get(theta_, rprime_)
  
    implicit none
    REAL(8), INTENT(IN) :: theta_(:), rprime_(:)
    real(8) :: thetaprime(N_PAR), rtemp, rtemp2, rtemp3, rtemp_wgt, rtemp_Pearson
    
    real(8) :: gmmse_WD(N_moments_all, N_PAR), gmmse_DWD(N_PAR, N_PAR)   
    integer(4) :: i, ii, imoments, j, k, jj, kk, kk10, kk2, kk20
    !this is for XI^2 statistics
    
    !====Now calculate the moments if parameters deviate one by one
    do i=1, N_PAR
       thetaprime(1:N_PAR) = theta_(1:N_PAR)
       if (gcL_partrans) then
          rtemp = real2ab_exp(gr_paradomain(i,1), gr_paradomain(i,2), theta_(i)) * rprime_(i)
          thetaprime(i) = ab2real_exp(gr_paradomain(i,1), gr_paradomain(i,2), rtemp)       
       else
          thetaprime(i) = theta_(i) * rprime_(i)
       endif
       
       rtemp = SMM_distance(thetaprime(1:N_PAR))       
       gmmse_simmoments%a(1:gmmse_simmoments%d1,1:MSM_Moments_Sim%d2,i)  &
             & = MSM_Moments_Sim%a(1:gmmse_simmoments%d1,1:MSM_Moments_Sim%d2)
    enddo   !do i
    
    ! moments:
    ! 1 lfpr_mean, 2 lnw_mean, 3 lnw_sd, 4 lnw_fe_mean, 
    ! 5 consumption  6 ss,  7 diff_lfpr_E2G, 8 diff_lfpr_G2B, 9 diff_lfpr_B2D 
    ! 1.1 E to U, 2.1 U to E, 3.1 wage depreciation after unemp
    
    !====calculate the variance-covariance matrix V. French&Jones(2011), supplement, page 12.
    !====Calculate the Jacobian matrix of population moment vector: #moment X #paras
          
    gmmse_jacobian%a = 0.0d0
    imoments = 0
    do ii = 1, N_moments_type
       kk10 = N_data_10
       if (ii == 5) then  !5 consumption, from 27
          kk = N_data_0_c 
       elseif (ii == 6) then  !6 ss
          kk = N_data_0_ss
          kk10 = N_data_10_ss
       elseif ((ii>=7) .AND. (ii<=10)) then  !#7 #8 #9 diff_lfpr
          if (IHasHealthAge>1 .AND. ((ii>=7) .AND. (ii<=9)) ) then
              kk = N_data_0_health
          elseif (Ipt==1 .and. (ii==7 .OR. ii==10)) then  !PT
              kk = N_data_0
          else
              cycle 
          endif          
       else
          kk = N_data_0  !22
       endif
       do j = kk, kk10
          imoments = imoments + 1
          do k = 1, N_PAR
             if (gcL_partrans) then 
                rtemp = real2ab_exp(gr_paradomain(k,1), gr_paradomain(k,2), theta_(k)) 
             else
                rtemp = theta_(k)
             endif
             gmmse_jacobian%a(imoments,k) = gmmse_simmoments%a(j,ii,k) - gmmse_simmoments%a(j,ii,gmmse_simmoments%d3)
             gmmse_jacobian%a(imoments,k) = gmmse_jacobian%a(imoments,k) / (rtemp*(rprime_(k)-1.0d0))  !numerical derivative
          enddo !k
       enddo !j
    enddo !ii
    !1.1 E to U, 2.1 U to E, 3.1 wage depreciation after unemp
    do j = 1, 3
       imoments = imoments + 1 
       do k = 1, N_PAR
          if (gcL_partrans) then 
             rtemp = real2ab_exp(gr_paradomain(k,1), gr_paradomain(k,2), theta_(k)) 
          else
             rtemp = theta_(k)
          endif
          gmmse_jacobian%a(imoments,k) = gmmse_simmoments%a(j,1,k) - gmmse_simmoments%a(j,1,gmmse_simmoments%d3)
          gmmse_jacobian%a(imoments,k) = gmmse_jacobian%a(imoments,k) / (rtemp*(rprime_(k)-1.0d0))  !numerical derivative
       enddo !k
    enddo !j       
             
    !write the Jacobian matrix
    open(15,file=trim(cwd)//'output/txt_par_jacobian'//trim(gmmse_filename)//'.txt')
    do i = 1, gmmse_jacobian%d1
       do j = 1, gmmse_jacobian%d2
          write(15,'(I4,A,I4,A,E20.10)')  i, ', ', j, ', ', gmmse_jacobian%a(i,j)
       enddo
    enddo
    close(15)
    
    if (gcI_robust == 21) then !no H depreciation at all
        ii = N_moments_all-1
        jj = N_PAR - 1
        do i = 2, N_PAR - 1
            gmmse_jacobian%a(:,i) = gmmse_jacobian%a(:,i+1)
        enddo
        
        gmmse_WD(1:(imoments-1),1:jj) = MATMUL(gmmse_W%a(1:ii,1:ii), gmmse_jacobian%a(1:ii,1:jj))    
        gmmse_DWD(1:jj,1:jj) = MATMUL(TRANSPOSE(gmmse_jacobian%a(1:ii,1:jj)), gmmse_WD(1:ii,1:jj))
        gmmse_DWD(1:jj,1:jj) = Matrix_Inverse(gmmse_DWD(1:jj,1:jj))
        gmmse_WD(1:imoments-1,1:jj) = MATMUL(gmmse_WD(1:ii,1:jj), gmmse_DWD(1:jj,1:jj)) !WD(D'WD)^{-1}    
        gmmse_DWD(1:jj,1:jj) = MATMUL(TRANSPOSE(gmmse_WD(1:ii,1:jj)), MATMUL(gmmse_S%a(1:ii,1:ii), gmmse_WD(1:ii,1:jj)))  !V 
        
        do i = N_PAR, 3, -1
            gmmse_DWD(i,i) = gmmse_DWD(i-1,i-1)
        enddo
        gmmse_DWD(2,2) = 0.0d0  !depreciation
        
    else
        gmmse_WD(1:imoments,:) = MATMUL(gmmse_W%a, gmmse_jacobian%a)    
        gmmse_DWD = MATMUL(TRANSPOSE(gmmse_jacobian%a), gmmse_WD)
        gmmse_DWD = Matrix_Inverse(gmmse_DWD)
        gmmse_WD(1:imoments,:) = MATMUL(gmmse_WD, gmmse_DWD) !WD(D'WD)^{-1}    
        gmmse_DWD = MATMUL(TRANSPOSE(gmmse_WD), MATMUL(gmmse_S%a, gmmse_WD))  !V       
    endif   !if (gcI_robust == 21)
    
    rtemp = dble(N_Data) / dble(N_Sim)  !\tau
    rtemp2 = (1.0d0+rtemp)/dble(N_Data)
    rtemp = sqrt(rtemp2)
    open(19,file=trim(cwd)//'output/txt_par_se'//trim(gmmse_filename)//'.txt')
    do i=1, N_PAR
        write(19,'(I4,A,E20.10,A,E20.10)')  i, ', ', gmmse_DWD(i,i)*rtemp2,  &
                              &  ', ', sqrt(abs(gmmse_DWD(i,i)))*rtemp
        !write(*,*)  '(var,se) for par # ',i, ' is (', gmmse_DWD(i,i)*rtemp2,  &
        !                       &  ', ', sqrt(abs(gmmse_DWD(i,i)))*rtemp, ')'
    enddo
    
    !!!!!!!!!!!!!!!! calculate Chi^2 statistics !!!!!!!!!!!!!!!
    ! First of all, get Matrix R
    if (gcI_robust == 21) then !no H depreciation at all
        gmmse_WD(1:ii,1:jj) = MATMUL(gmmse_W2%a(1:ii,1:ii), gmmse_jacobian%a(1:ii,1:jj))
        gmmse_DWD(1:jj,1:jj) = MATMUL(TRANSPOSE(gmmse_jacobian%a(1:ii,1:jj)), gmmse_WD(1:ii,1:jj))
        gmmse_DWD(1:jj,1:jj) = Matrix_Inverse(gmmse_DWD(1:jj,1:jj))
        gmmse_XI2_PR%a(1:ii,1:ii) = - MATMUL( MATMUL(gmmse_jacobian%a(1:ii,1:jj), gmmse_DWD(1:jj,1:jj)),  &
                                        &     TRANSPOSE(gmmse_WD(1:ii,1:jj)) )
        do i = 1, ii
            gmmse_XI2_PR%a(i,i) = 1.0d0 + gmmse_XI2_PR%a(i,i)
        enddo
        gmmse_XI2_PR%a(1:ii,1:ii) = MATMUL(gmmse_XI2_PR%a(1:ii,1:ii), MATMUL(gmmse_S%a(1:ii,1:ii),  &
                                        &  gmmse_XI2_PR%a(1:ii,1:ii))) 
        gmmse_XI2_PR%a(1:ii,1:ii) = Matrix_Inverse(gmmse_XI2_PR%a(1:ii,1:ii))
    else    
        gmmse_WD = MATMUL(gmmse_W2%a, gmmse_jacobian%a)
        gmmse_DWD = MATMUL(TRANSPOSE(gmmse_jacobian%a), gmmse_WD)
        gmmse_DWD = Matrix_Inverse(gmmse_DWD)
        gmmse_XI2_PR%a = - MATMUL( MATMUL(gmmse_jacobian%a, gmmse_DWD), TRANSPOSE(gmmse_WD) )
        do i = 1, gmmse_XI2_PR%d1
            gmmse_XI2_PR%a(i,i) = 1.0d0 + gmmse_XI2_PR%a(i,i)
        enddo
        gmmse_XI2_PR%a = MATMUL(gmmse_XI2_PR%a, MATMUL(gmmse_S%a, gmmse_XI2_PR%a)) 
        gmmse_XI2_PR%a = Matrix_Inverse(gmmse_XI2_PR%a)
    endif
       
    
    !Now calculate the Chi^2 statistic and Pearson's test
    rtemp = 0.0d0
    rtemp_Pearson = 0.0d0  !Pearson's test
    imoments = 0  !second dimension
    do i = 1, N_moments_type
       kk10 = N_data_10 
       if (i == 5) then  !5 consumption, from 27
          kk = N_data_0_c 
       elseif (i == 6) then !6 ss
          kk = N_data_0_ss
          kk10 = N_data_10_ss
       elseif ((i>=7) .AND. (i<=10) ) then  !#7 #8 #9 diff_lfpr
          if (IHasHealthAge>1 .AND. ((i>=7) .AND. (i<=9)) ) then
              kk = N_data_0_health
          elseif (Ipt==1 .and. (i==7 .OR. i==10)) then  !PT
              kk = N_data_0
          else
              cycle 
          endif          
       else
          kk = N_data_0  !22
       endif
       
       do ii = kk, kk10
          imoments = imoments + 1
          
          !Pearson's test
          rtemp_Pearson = rtemp_Pearson + (gmmse_simmoments%a(ii,i,gmmse_simmoments%d3)-MSM_Moments_Data%a(ii,i))**2 &
                                       &  / MSM_Moments_Data%a(ii,i)
          
          !Chi^2 statistic
          k = 0 !first dimension
          rtemp3 = 0.0d0
          do j = 1, N_moments_type         
             kk20 = N_data_10 
             if (j == 5) then  !5 consumption, from 27
                kk2 = N_data_0_c 
             elseif (j == 6) then  !6 ss
                kk2 = N_data_0_ss
                kk20 = N_data_10_ss
             elseif ((j>=7) .AND. (j<=10) ) then  !#7 #8 #9 diff_lfpr
                 if (IHasHealthAge>1 .AND. ((j>=7) .AND. (j<=9)) ) then
                     kk2 = N_data_0_health
                 elseif (Ipt==1 .and. (j==7 .OR. j==10)) then   !PT
                     kk2 = N_data_0
                 else
                     cycle 
                 endif
             else
                kk2 = N_data_0  !22
             endif
              
             do jj = kk2, kk20
                k = k + 1
                !rtemp_wgt = gmmse_XI2_PR%a(k,imoments)
                rtemp_wgt = gmmse_W2%a(k,imoments)
                rtemp3 = rtemp3 + (gmmse_simmoments%a(jj,j,gmmse_simmoments%d3)-MSM_Moments_Data%a(jj,j))*rtemp_wgt
             enddo !jj
          enddo !j
          !1.1 E to U, 2.1 U to E, 3.1 wage depreciation after unemp
          do jj = 1, 3
              k = k + 1
              
              if (gcI_robust==21 .AND. jj==3) cycle     !no H depreciation at all
              
              rtemp_wgt = gmmse_W2%a(k,imoments)
              rtemp3 = rtemp3 + (gmmse_simmoments%a(jj,1,gmmse_simmoments%d3)-MSM_Moments_Data%a(jj,1))*rtemp_wgt
          enddo !j
          
          rtemp = rtemp + rtemp3 * (gmmse_simmoments%a(ii,i,gmmse_simmoments%d3)-MSM_Moments_Data%a(ii,i))
       enddo !ii
    enddo ! i
    
    !1.1 E to U and 2.1 U to E, 3.1 wage depreciation after unemp
    i = 1
    do ii = 1, 3
       imoments = imoments + 1
       
       if (gcI_robust==21 .AND. ii==3) cycle     !no H depreciation at all 
       
       !Pearson's test
       rtemp_Pearson = rtemp_Pearson + (gmmse_simmoments%a(ii,i,gmmse_simmoments%d3)-MSM_Moments_Data%a(ii,i))**2 &
                                    &  / MSM_Moments_Data%a(ii,i)
       !Chi^2 statistic
       k = 0 !first dimension
       rtemp3 = 0.0d0
       do j = 1, N_moments_type
          kk20 = N_data_10 
          if (j == 5) then  !5 consumption, from 27
             kk2 = N_data_0_c 
          elseif (j == 6) then  !6 ss
             kk2 = N_data_0_ss
             kk20 = N_data_10_ss
          elseif ((j>=7) .AND. (j<=10) ) then  !#7 #8 #9 diff_lfpr
              if (IHasHealthAge>1 .AND. ((j>=7) .AND. (j<=9)) ) then
                  kk2 = N_data_0_health
              elseif (Ipt==1 .and. (j==7 .OR. j==10)) then   !PT
                  kk2 = N_data_0
              else
                  cycle 
              endif
          else
             kk2 = N_data_0  !22
          endif
              
          do jj = kk2, kk20
             k = k + 1
             !rtemp_wgt = gmmse_XI2_PR%a(k,imoments)
             rtemp_wgt = gmmse_W2%a(k,imoments)
             rtemp3 = rtemp3 + (gmmse_simmoments%a(jj,j,gmmse_simmoments%d3)-MSM_Moments_Data%a(jj,j))*rtemp_wgt
          enddo !jj
       enddo !j
       !1.1 E to U and 2.1 U to E, 3.1 wage depreciation after unemp
       do jj = 1, 3
          k = k + 1
          
          if (gcI_robust==21 .AND. jj==3) cycle     !no H depreciation at all 
          
          rtemp_wgt = gmmse_W2%a(k,imoments)
          rtemp3 = rtemp3 + (gmmse_simmoments%a(jj,1,gmmse_simmoments%d3)-MSM_Moments_Data%a(jj,1))*rtemp_wgt
       enddo !j
       
       rtemp = rtemp + rtemp3 * (gmmse_simmoments%a(ii,i,gmmse_simmoments%d3)-MSM_Moments_Data%a(ii,i))
    enddo !ii

    
    rtemp = rtemp / rtemp2
    write(19,'(A,I5,A,I5,A,I5)')  'M=', imoments, ', n=', N_PAR, ', M-n=', imoments-N_PAR
    write(19,'(A,E50.10)')  'CHI^2=', rtemp
    write(19,'(A,E50.10)')  'Pearson=', rtemp_Pearson
    close(19)
    !!!!!!!!!!!!!!!!ENDOF calculate XI^2 statistics
    
    return

  end subroutine sub_GMM_SE_get
  
  ! varying one parameter and see how the results change
  subroutine sub_GMM_onepar(theta_, ipar_)  !called in the master
    implicit none
    REAL(8), INTENT(IN) :: theta_(:)
    integer, INTENT(IN) :: ipar_
    integer(4), parameter :: N_onepar = 21  ! number of varying the parameter value.
    real(8) :: thetaprime(N_PAR), r0, rs, rVtemp(N_onepar), rtemp, r1, r2, r3, r4
    integer(4) :: i, j, i1, i2
  
    gmmse_simmoments%d1 = N_T  !max(N_data_10, N_data_10_ss)
    gmmse_simmoments%d2 = MSM_Moments_Data%d2  !6 (no health) or 9 (with health) or 7 (pt only) or 10 (pt and health)
    gmmse_simmoments%d3 = N_onepar  !the first one is the original one
    allocate(gmmse_simmoments%a(gmmse_simmoments%d1, gmmse_simmoments%d2, gmmse_simmoments%d3))
    gmmse_simmoments%a = 0.0d0
    
    ! moments:
    ! 1 lfpr_mean, 2 lnw_mean, 3 lnw_sd, 4 lnw_fe_mean, 
    ! 5 consumption  6 ss,  7 diff_lfpr_E2G, 8 diff_lfpr_G2B, 9 diff_lfpr_B2D, 7/10 pt
    ! 1.1 E to U, 2.1 U to E, 3.1 wage depreciation after unemp
        
    j = 4 !-0.77 -0.05*(1-4)
    
    !keep original moments in j+1. This has to run immediately after the output
    gmmse_simmoments%a(1:gmmse_simmoments%d1, 1:MSM_Moments_Sim%d2, j+1)  &
          & = MSM_Moments_Sim%a(1:gmmse_simmoments%d1,1:MSM_Moments_Sim%d2)
    
    thetaprime = theta_   !original parameters
    
    r0 = real2ab_exp(gr_paradomain(ipar_,1), gr_paradomain(ipar_,2), theta_(ipar_))  !-0.77 in the baseline model
    !toward lower end
    rs = 0.05d0
    do i = 1, j
        rVtemp(i) = r0 - rs*dble(j+1 - i) 
        thetaprime(ipar_) = ab2real_exp(gr_paradomain(ipar_,1), gr_paradomain(ipar_,2), rVtemp(i)) 
        rtemp = SMM_distance(thetaprime(1:N_PAR))       
        gmmse_simmoments%a(1:gmmse_simmoments%d1,1:MSM_Moments_Sim%d2,i)  &
             & = MSM_Moments_Sim%a(1:gmmse_simmoments%d1,1:MSM_Moments_Sim%d2)
    enddo
    
    rVtemp(j+1) = r0
    !toward upper end, same step of 0.05
    !rs = 0.05d0
    do i = j+2, N_onepar
        rVtemp(i) = min(0.0d0, r0 + rs*dble(i - (j+1))) 
        thetaprime(ipar_) = ab2real_exp(gr_paradomain(ipar_,1), gr_paradomain(ipar_,2), rVtemp(i))
        rtemp = SMM_distance(thetaprime(1:N_PAR))       
        gmmse_simmoments%a(1:gmmse_simmoments%d1,1:MSM_Moments_Sim%d2,i)  &
             & = MSM_Moments_Sim%a(1:gmmse_simmoments%d1,1:MSM_Moments_Sim%d2)
    enddo
    
    open(15,file=trim(cwd)//'output/txt_deviation_onepar.txt')
    if (Ipt==1) then   !PT
        i1 = N_data_0_pt
        i2 = max(N_data_10, N_data_10_ss, N_data_10_pt)
    else
        i1 = N_data_0
        i2 = max(N_data_10, N_data_10_ss)
    endif
        
    do i = 1, N_onepar
        do git = i1,  i2
            if (IHasHealthAge > 0) then
                r1 = gmmse_simmoments%a(git,7, i)  !diff_lfpr_E2G
                r2 = gmmse_simmoments%a(git,8, i)  !diff_lfpr_G2B
                r3 = gmmse_simmoments%a(git,9, i)  !diff_lfpr_B2D
            else
                r1 = -999.0d0
                r2 = -999.0d0
                r3 = -999.0d0
            endif
             
            if (Ipt == 1) then !PT
                if (IHasHealthAge > 0) then
                    r4 = gmmse_simmoments%a(git,10, i)  !pt
                else
                    r4 = gmmse_simmoments%a(git,7, i)  !pt
                endif
            else
                r4 = -999.0d0
            endif
           
            write(15,'(I8,A,F12.4,A,I4,A,F12.4,A,F12.4,A,F12.4,A,F12.4,A,F12.4, &
                   & A,F12.4,A,F12.4,A,F12.4,A,F12.4,A,F12.4,A,F12.4,A,F12.4,A,F12.4)')  &
                   & i, ',', rVtemp(i), ',', git,         &
                   & ',', gmmse_simmoments%a(git,1, i),   &
                   & ',', gmmse_simmoments%a(git,2, i),   &
                   & ',', gmmse_simmoments%a(git,3, i),   &
                   & ',', gmmse_simmoments%a(git,4, i),   &
                   & ',', gmmse_simmoments%a(1,1, i),     &   !E to U
                   & ',', gmmse_simmoments%a(2,1, i),     &   !U to E
                   & ',', gmmse_simmoments%a(git,5, i),   &   !consumption
                   & ',', gmmse_simmoments%a(git,6, i),   &   !#6 SS
                   & ',', r1, ',', r2, ',', r3, ',', r4,   &   !diff_lfpr_E2G, diff_lfpr_G2B, diff_lfpr_B2D, lfpr_pt
                   & ',', gmmse_simmoments%a(3,1, i)          !wage depreciation after unemployment spell

        enddo  !do git
    enddo ! do i    
          
    close(15)
    
    !deallocate variables
    deallocate(gmmse_simmoments%a)
    
  end subroutine sub_GMM_onepar

END MODULE GMM_SE
